The processing pipeline is build upon have a distinct raster for each Time series interpolation case, listing this into a single object and then using lapply to perform functions across each TSI case.
# base case
link <- "/data/Dagobah/greengrass/schnesha/03_RBF_interpolation_LND_1982-2021/base_case/mosaic/1984-2021_001-365_HL_TSA_LNDLG_NDV_TSS.vrt"
tsi_base_raw <- rast(link)
# case day_16_sigma_8_16_32
link <- "/data/Dagobah/greengrass/schnesha/03_RBF_interpolation_LND_1982-2021/day_16_sigma_8_16_32/mosaic/1984-2021_001-365_HL_TSA_LNDLG_NDV_TSI.vrt"
tsi_day_16_sigma_8_16_32_raw <- rast(link)
# case day_16_sigma_8_8_8_16_16_16_64
link <- "/data/Dagobah/greengrass/schnesha/03_RBF_interpolation_LND_1982-2021/day_16_sigma_8_8_8_16_16_16_64/mosaic/1984-2021_001-365_HL_TSA_LNDLG_NDV_TSI.vrt"
tsi_day_16_sigma_8_8_8_16_16_16_64_raw <- rast(link)
# case day_16_sigma_8_8_8_16_16_16_32_64
link <- "/data/Dagobah/greengrass/schnesha/03_RBF_interpolation_LND_1982-2021/day_16_sigma_8_8_8_16_16_16_32_64/mosaic/1984-2021_001-365_HL_TSA_LNDLG_NDV_TSI.vrt"
tsi_day_16_sigma_8_8_8_16_16_16_32_64_raw <- rast(link)
# All cases are combined into a single list to facilitate using lapply based processing
tsi_list <- list(tsi_base_raw, tsi_day_16_sigma_8_16_32_raw, tsi_day_16_sigma_8_8_8_16_16_16_64_raw, tsi_day_16_sigma_8_8_8_16_16_16_32_64_raw)
cases <- c("base_case", "day_16_sigma_8_16_32", "tsi_day_16_sigma_8_8_8_16_16_16_64", "tsi_day_16_sigma_8_8_8_16_16_16_32_64")
names(tsi_list) <- cases
tsi_obs_length <- list()
lapply(seq_along(tsi_list),
function(i)
tsi_obs_length[i] <<- length(names(tsi_list[[i]]))
)link <- "/data/Dagobah/greengrass/schnesha/03_RBF_interpolation_LND_1982-2021/DE_shp/vg2500_bld_EPSG3035.shp"
DE <- read_sf(link)create stack for direct location comparison unique index for each TSI case as number of layers differ between them
base_case_obs_count <- readRDS(file = "/data/Dagobah/greengrass/schnesha/03_RBF_interpolation_LND_1982-2021/R_objects/base_case_obs_count.rds")
tsi_stack <- c(tsi_list[[1]]$`20000514_LND05`,
tsi_list[[2]]$`20000509`,
tsi_list[[3]]$`20000509`,
tsi_list[[4]]$`20000509`)
plot(tsi_list[[1]]$`20160807_LND08`)
plot(vect(DE), add=TRUE)#de <- st_read("/data/Dagobah/dc/deu/vector/DEU.gpkg")
#plot(tsi_list[[2]][[800:808]])
for (i in seq_along(tsi_list)) {
plot(tsi_stack[[i]],
legend=T, col = viridis(n=100,option="D"),
range=c(0,10000),
ext = ext(DE),
plg=list(title= paste(names(tsi_list[i]), "\n", names(tsi_stack[[i]]))
))
plot(vect(DE), add=TRUE)
}can be enables if desired, but is computationally quite demanding
val = as.numeric(c(0:10000))
pal = colorNumeric(c("viridis"), val,
na.color = "grey")
#pal <- colorNumeric(c("#0C2C84", "#41B6C4", "#FFFFCC"), values(tsi_list[[4]][[675]]), na.color = "grey")
leaflet() |>
addTiles() |>
addRasterImage(tsi_list[[4]][[675]], group = "show_layer",
maxBytes = (20 * 1024 * 1024)) |>
addLegend(pal = pal,
values = values(tsi_list[[4]][[600:600]]),
title = "Reflectance") |>
addLayersControl(
overlayGroups = c("show_layer"),
options = layersControlOptions(collapsed = FALSE)
)## Warning in pal(c(r[1], cuts, r[2])): Some values were outside the color scale
## and will be treated as NA